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Abstract 

We introduce unconditionally stable finite element approximations for a phase 
field model for solidification, which take highly anisotropic surface energy and ki- 
netic effects into account. We hence approximate Stefan problems with anisotropic 
Gibbs-Thomson law with kinetic undercooling, and quasi-static variants thereof. 
The phase field model is given by 

$w t + \ g{<p) (ft = V . (b(ip) V w) , 

c* I g(<p) w = e £ /i(V <p) <pt - e V . A'(V if) + e' 1 tf'(p) 

subject to initial and boundary conditions for the phase variable (p and the temper- 
ature approximation w. Here e > is the interfacial parameter, ^ is a double well 
potential, = \/2 W(s) ds, g is a shape function and A(V <p) = \ |7(Vc^)| 2 , 
where 7 is the anisotropic density function. Moreover, # > 0, A>0, a > 0, a > 
and p > are physical parameters from the Stefan problem, while b and \x are 
coefficient functions which also relate to the sharp interface problem. 

On introducing the novel fully practical finite element approximations for the 
anisotropic phase field model, we prove their stability and demonstrate their appli- 
cability with some numerical results. 
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1 Introduction 

Phase field models are a successful approach for interface evolution in cases where inter- 
facial energy is important, and many numerical approaches for the underlying equations 
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have been studied in the literature. However, in situations where anisotropy is incorpo- 
rated only very few results related to the numerical analysis of approximations to the 
phase field system have appeared in the literature. The reason for this is that the un- 
derlying equations involve highly nonlinear parabolic partial differential equations. Since 
phase field models describe very unstable solidification phenomena, it seems to be very 
important to use stable approximation schemes which do not trigger additional instabil- 
ities resulting from discretization errors. In this context we would like to mention that 
there exist many computations on anisotropic solidification, with the help of phase field 
equations, showing pattern formation which is driven by the discretization rather than 
by the underlying partial differential equations. The goal of this paper is to introduce 
and analyze a new stable finite element approximation for the anisotropic phase field sys- 
tem. The approach is based on earlier work for the Allen-Cahn and the Cahn-Hilliard 
equations, see Barrett et al. (2012c), and on ideas on how to handle the anisotropy that 



have been used earlier for sharp interface models by the same authors, see |Barrett et al. 
Q2008a|bD . To our knowledge, the introduced finite element approximation is the first 
unconditionally stable approximation of a phase field model for anisotropic solidification 
in the literature. 

As the phase field model and its quasi-stationary variant, the viscous Cahn-Hilliard 
equation, converge to sharp interface models for solidification in the asymptotic limit as 
the interfacial thickness tends to zero we first introduce the sharp interface model. Let 
Y{t) C M. d , d = 2,3, denote the interface between a solid and liquid phase, say, or a solid 
phase and a gas phase. Then the surface energy of T(t) is defined as 



7(n) ds. 



r(t) 



where n denotes the unit normal of T(t), and where the anisotropic density function 
7 : R d — > ]R> with 7 e C 2 (R d \ {0}) n C(M. d ) is assumed to be absolutely homogeneous 
of degree one, i.e. 

7 (Ap) = |A|7(p) VpeR^VXeR = 

with 7' denoting the gradient of 7. 



Y(p).p = 7 (p) VpeM d \{0}, (1.2) 



Relevant for our considerations is the first variation, — k 7 , of (1.1), which can be 
computed as 

k 7 := -V s . 7 ; (n) ; 



where V s . is the tangential divergence of T, see e.g. Cahn and Hoffman (1974); Barrett 
et al. (2008b, 2010b). Note that k 7 reduces to the sum of the principal curvatures of T in 



the isotropic case, i.e. when 7 satisfies 

lip) = \p\ 



;i.3) 



Then the full Stefan problem that we want to consider in this paper is given as follows, 
where Q C R d is a given fixed domain with boundary dQ and outer normal v. 
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Find u : O X [0, T] — > R and the interface (T(t)) t ^[o t T] such that for all t G (0, T] the 
following conditions hold: 

tfu t -lC-Au = inft_(t), $u t -]C + Au = 
du 



dn 



-XV 



T(t) 

/3(n) 
du 
du 

r(o) = r , 



a — au 



on c^fi, u = ud 



in Q + (t), 


(1.4a) 


on r(t), 


(1.4b) 


on r(t), 


(1.4c) 


on c^ft , 


(1.4d) 


in ft . 


(1.4e) 



In the above u denotes the deviation from the melting temperature Tm, i.e. Tm is the 
melting temperature for a planar interface. In addition, ft_(i) is the solid region, with 
boundary T(t) = <9ft_(t), so that the liquid region is given by ft+(t) := ft \ ft_(t). Here 
we assume that the solid region ft_(£) has no intersection with the external boundary 
dfl, but more general situations can also be considered, as will be outlined in Section [3] 
below. Moreover, here and throughout this paper, for a quantity v defined on Q, we use 
the shorthand notations t>_ := v \q_ and v + := v \q + . The parameters $ > 0, A > 0, 
p>0, a>0, a > are assumed to be constant, while K,± > are assumed to be 
constant in each phase. The mobility coefficient (3 : M. d — > M> is assumed to satisfy 
f3(p) > for all p ^ and to be positively homogeneous of degree one. In addition 
[JC g]r(t)(z) := (/C + ^ - K- ^){z) for all z E T{t), and V is the velocity of T(t) in the 
direction of its normal n, which from now on we assume is pointing into Q + (t). Finally, 
dQ = djyQ U dpQ with d^Q fl dpQ = 0, up : c^fi — )■ R is the applied supercooling at the 
boundary, and r C Q and Mq : ^ ~^ R are given initial data. 



The model (1.4a-e) can be derived for example within the theory of rational ther- 
modynamics and we refer to Gurtin (1988) for details. We remark that a der ivation 
from thermodynamics would lead to the identity a — We note that (1.4b) is the 
well-known Stefan condition, while (1.4c) is the Gibbs-Tliomson condition, with kinetic 



undercooling if p > 0. The case i9>0,p>0,a:>0 leads to the Stefan problem with the 
Gibbs-Thomson law and kinetic undercooling. In some models in the literature, see e.g. 
Luckhaus (1990), the kinetic undercooling is set to zero, i.e. 



p = 0. Setting ■d = p = 
but keeping a > leads to the Mullins-Sekerka problem with the Gibbs-Thomson law, 

sec 



Mullins and Sekerka (1963). 



For later reference, we introduce the function spaces 

S ■= {rj G : r] = on d D n} and S D := {r] G H 1 ^) : rj = u D on d D Vl} 

where we assume for simplicity of the presentation from now on that 

either (i) dVL = d^Vl , (ii) dVl = , 

or (iii) ft = (-H, H) d } d D n = [-H, H}*- 1 x {H}, H>0; 



(1.5) 
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and, in the cases (1.5) (i) and (iii), that ud G Hs(djjQ). For notational convenience, we 



define ud '■= in the case (1.5) (ii 



We recall from Barrett et al. (2010b) that, on assuming that ujj is constant, for a 



solution u and V to (1.4a-e) it can be shown that the following formal energy equality 
holds 



dt V 2 



\u - u D \ + 



A a 



7(n) ds - Xu D \n+(t)\ + (K V«, Vu 



r(t) 



+ 



Ap 



V 2 

\\t) /5(n) 



ds = 0, (1.6) 



where 



denotes the L 2 -inner product over Q, with the corresponding norm given by 



| ■ |o, and where := J n+ ^ 1 dx. 

In Section [2] we will precisely state a phase field model which approximates the free 
boundary problem ( 1.4a{ -e). We only mention here that the phase field method is based 
on the idea of a diffuse interface, which hence has a positive thickness. Let us briefly 
discuss some relevant literature. For solidification the phase field method was originally 
proposed by Langer (1986) as a model for solidification of a pure substance. It was 



Kobayashi (1993) who first was able to simulate complicated dendritic patterns which 



resemble those appearing during solidification. Since then an enormous effort has gone 
into numerically studying phase field models. We refer only to Elliott and Gardiner 



(1996) 



(2002) 



Karma and Rappel (1996, 1998) and to the reviews Boettinger et al. (2002); Chen 



McFadden (2002); Singer-Loginova and Singer (2008). 



(1.4a 



A phase field model, and its numerical approximation, for the sharp interface problem 
with $ = p = and fC + = XL has been considered in the recent paper Barrett 



et al. 



(2012c). In particular, the authors were able to present unconditionally stable finite 
where the treatment of the anisotropy does not lead to new 

It is one of the aims of the present 



element approximations 

nonlinearities compared to the isotropic situation. 



article to extend the discretizations in Barrett et al. (2012c) to the more general problem 
(1.4a-e), i.e. in particular to the case $ > 0, and p > 0, and to a wider class of anisotropies 



than considered in Barrett et al. (2012c). The new anisotropies considered in the present 
article will lead to more nonlinear schemes, however. 

The remainder of the paper is organized as follows. In Section [2] we state the two phase 
field models for the approximation of the sharp interface problem (1.4a-e) that we want 



to consider in this paper. In Section [3] we introduce our finite element approximations for 
these problems, and we prove stability results for these approximations. Solution methods 
for the discrete equations are shortly reviewed in Section |4j In addition, we present several 
numerical experiments in Section [5} 
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2 Phase field models and anisotropics 



Phase field models are a computational tool to compute approximations for sharp interface 
evolutions such as (1.4a-e), without having to capture the sharp interface T(t) directly. 
On introducing a phase field tp : Q x (0, T) — > R, where the sets fi±(£) := {x G Q : 
±<p(x,t) > 0} are approximations to fl±(t), a system of partial differential equations for 
(f can be derived so that the zero level sets of ip formally approximate the interface T(t), 
satisfying e.g. (1.4a-e), in a well defined limit. For more details on phase field methods 



and other approaches to the approximation of the evolution of interfaces we refer to the 



review article Deckelnick et al. (2005) and the references therein. 



On introducing the small interfacial parameter e > 0, it can be shown that 

— S^((p) « / 7(n) ds, 

for e sufficiently small, where 

£ 7 (^) : = ! fMWOf + e-^foOdx with c * : =/ y^{s) ds . (2.1) 
Jn J-i 

Here ¥ : R — > [0, oo] is a double well potential, which for simplicity we assume to be 
symmetric and to have its global minima at ±1. The canonical example is 



¥00 := \ (s 2 - l) 2 
Another possibility is to choose 



*'(s) = s 3 - s and Cvj> = | 2§ . 



oo 



l s l < 1 

\s\ > 1 . 



7T . 

2 ' 



(2.2) 



(2.3) 



see e.g. |Blowey and Elliot^ ( |1992| ); |EUiott and Gardiner| p996[ ); |Elliott| p997| ). Clearly 
the obstacle potential (2.3), which forces (p to stay within the interval [—1,1], is not 
differentiable at ±1. Hence, whenever we write \l/'(s) in the case (2.3) in this paper, we 
mean that the expression holds only for \s\ < 1, and that in general a variational inequality 
needs to be employed. While it can be shown that the asymptotic interface thickness in 



phase field models with (2.1) for the isotropic surface energy (1.3) is proportional to e 



for anisotropic energy densities the asymptotic interface thickness is no longer uniform, 



but now also depends on 7 and on V (p, see e.g. Bellettini and Paolini (1996); Elliott and 



Schatzle (1996); Wheeler and McFadden (1996) 



We remark that other, non-classical, phase field models are based on the energy 

/ \V<p\~ 1 'y{V<p) (I \Vip\ 2 + e- 1 ^(i P )) dx 
Jci 



(2.4) 



for e.g. the smooth double-well potential (2.2), see Torabi et al. (2009). The energy (2.4) 



has the advantage that the asymptotic interface thickness is now only determined by e 
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(independently of 7 and the orientation of the interface), whereas the disadvantage is 
that the resultant partial differential equations become more nonlinear and are singular 
at Vy? = 0. We note that higher order regularizations of the energies (2.1) and (2.4) 



in the case of a non-convex anisotropy density function 7, which lead to sixth order 



Cahn-Hilliard type equations, have been considered in e.g. Li et al. (2009). 



We are not aware of any numerical analysis for discretizations of anisotropic phase 



field models for (1.4a-e) involving either (2.1) or (2.4) 



We now state the two phase field models that we are going to consider in this paper. 
To this end, for p G IR d , let 



A{p) = \\i{p)\' 



and define 



A'(p) 
l{p) 



7(p)V(p) 
p = 



a p 



p^O, 

0. 



(2-5) 



(2-6) 



where /i G R is a constant satisfying min p ^ 



7 0) 

/3(P) 



< U < 



\x s maXp^o 



7(p) 



2.1 Viscous Cahn— Hilliard equation 



in 



A phase field model for (1.4a-e) with •& = p = has been recently studied by the authors 
Barrett et al. (2012c). The case $ = and p > gives rise to the following viscous 



Cahn-Hilliard equation for the anisotropic Ginzburg-Landau energy (2.1), where if is a 
phase field approximation to the (rescaled) temperature u: 



— w 
a 

dv 
w 
dw 
dv 



b(<p) 



V.(%)Vw) 

e £ //(V <p) (p t - e V . A'(V (p) + e~ l V'(<p) 
a 

0. 

u D 

0, 
^0 



in 


n T := 


:flx(o,r), 


(2.7a) 


in 


fir, 




(2.7b) 


on 


dtt x 


(0,T), 


(2.7c) 


on 


o D n 


x (0,T), 


(2.7d) 


on 




x (0,T), 


(2.7e) 


in 






(2.7f) 



where 



b(s) 



|(l + s)/C + + i(l-s) /C_ 



(2-8) 



With the help of formal asymptotics, see e.g. McFadden et al. (1993); Wheeler and Mc 



Fadden (1996); Bellettini and Paolini (1996), it can be shown that the sharp interface 



limit of (2.7a-f), i.e. the limit as e — > 0, is given by the quasi-static Stefan problem (or 
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Mullins-Sekerka problem) (1.4a-e) with •& = 0, and with u denoting the sharp interface 
limit of w. 



We remark that the phase field analogue of the sharp interface energy identity (1.6) 
in the case d = is given by the formal energy bound 



d ( \ a 1 
dt 



— —E^-IXud [ ipdx]+(b(if)Vw,Vw)+E^ — (fi(Vip),(ip t ) 2 )<0 
a c<j, Jn / a 

(2.9) 

for the phase field model (2.7a-f) with the potential (2.3). For smooth potentials such as 
(2.2) the energy law (2.9) holds with equality. 



2.2 Heat equation coupled to Allen— Cahn 



The second phase field model is based on the work in Elliott and Gardiner (1996), see also 



Kobayashi (1993); Karma and Rappel (1996) for other related approaches, and allows the 



sharp interface limit (1.4a-e) with $ > 0. It consists of a heat equation for the phase field 
temperature approximation w coupled to an Allen-Calm phase field equation for ip. In 
particular, we have the modified heat equation 

(2.10a) 

T) , (2.10b) 
T), (2.10c) 
(2.10d) 

where b is defined in (2.8), and where the function g G C 1 (1R) is such that 

g(s)>0 VsG[-l,l], J Q{y)&y = l and P(s) := J ' e {y) dy . 

We note that P, which is a monotonically increasing function over the interval [—1,1] 
with P(— 1) = and P(l) = 1, is often called the interpolation function. In this paper, 
we follow the convention from Elliott and Gardiner (1996), where g = P' is called the 



$w t + \ g(<f) <f t 


= V.(%)Vw) 


in Qt 5 


w 


= u D 


on do®- x 




= 


on x 




= i9 wq 


in Q , 



shape function. More details on interpolation functions P, respectively shape functions 



g, can be found in e.g. Wang et al. (1993); Garcke and Stinner (2006); Eck et al. (2006) 



Caginalp et al. (2008). In particular, if one also assumes symmetry, i.e. 

1,1], 



g{s) 



g(-s) V s G 



then a faster convergence of the phase field model to the sharp interface limit, as e — > 0, 
can be shown on prescribing suitable first order corrections in e for the remaining phase 



field parameters; see Karma and Rappel (1996, 1998); Eck et al. (2006); Garcke and 



Stinner (2006) for details. Possible choices of g that will be considered in this paper are 

2 



i 

2 ' 



ill gis 



1(1-*) 



(m) g{s) 



15 
16 



[s 2 -D 



(2.H) 



The heat equation (2.10a-d) is coupled to the following modified Allen-Calm equation: 

a , , p 



c>i> - gUp) w = e- /i(V <p) (pt - e V . A'(V <p) + £ V'(<p) 
a a 

av 



in Qi 



(2.12a) 



on dVl x (0,T) , (2.12b) 
in ft. (2.12c) 



We remark that the phase field analogue of the sharp interface energy identity ( 1.6[ ) is 
given by the formal energy bound 



d (•& . l2 \a 1 

at \ 2 a 



E^-Xud j P(<f) dx^j + (%) Vw, Vw) 



+ £ ^-(Mv v ),(^ 2 )<o 



(2.13) 



for the phase field model (2.10a-d), (2.12a-c) with the potential (2.3). For smooth poten- 



tials such as (2.2) the energy law ( |2.13 ) holds with equality. We remark that the energy 
decay in (2.13) for the phase field model (2.10a-d), (2.12a-c) means that the model can 



be said to be thermodynamically consistent. For more details on thermodynamically 



consistent phase field models we refer to e.g. Penrose and Fife (1990); Wang et al. (1993). 



Remark. 2.1. We remark that in the special case $ = ; and if we choose (2.11 )(i), 
then clearly ( 2.10af -d), (2.12a-c) collapses to the system (2.7a-f). Similarly, the energy 
law (2.13) in this case collapses to (2.9). Hence from now on in this paper, we will only 



consider the more general model (2.10a-d), (2.12a-c). Finally we note that the phase field 



model (2.7a-f) in the case p = was recently considered in Barrett et al. (2012c). 



We observe that for e small, on recalling that the thickness of the interfacial region 
goes to zero as e — > 0, it holds that 



p(p) dx « p(i) \nut)\ + p(-i) \n e _(t)\ = \nut)\ 



(2.14) 



which is a consequence of the fact that P(<p) approximates the characteristic function 



of the liquid phase ft + (i). It is clear from (2.13) and (2.14) that for negative values of 
Ud, V 9 is encouraged to take on negative values, so that the approximate liquid region 
ft+(i) shrinks, whereas positive values of ud encourage ip to take on positive values, so 
that the liquid region grows. Of course, this is simply the phase field analogue of the 



sharp interface behaviour induced by (1.6). A side effect of the interpolation function P 



in (2.13), however, is that the function 



G(s) = a (acq,e) ^(s)—u D P(s 



(2.15) 



need no longer have local minima at s — ±1. This can result, for example, in undesired, 
artificial boundary layers for strong supercoolings, i.e. when —ud is large; see also Re- 
marks 3^ and 3.9| below. For smooth potentials ^, sufficient conditions for s = ±1 to 
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be local minimum points of G(s) are g{+l) = f?'(±l) = 0, which is evidently satisfied by 



(2.11 ) (iii) - In fact, in applications phase field models for solidification almost exclusively 



use the quartic potential (2.2) together with this shape function; see e.g. Boettinger et al 



(2002 


); 


Chen 


(2002 


); 


McFadden 


(2002) 



For the obstacle potential (2.3) the situation is similar, although there is more flex- 
ibility in the possible choices of g. In particular, here a sufficient condition for G(s) to 
have local minima at s = ±1 is given by 



a(acys) 1 ±udq(±1)>0. 



(2.16) 



Clearly, (2.16) is always satisfied for (2.11 ) (iii) 

q0) 



while for ud < it is sufficient to require 
0, e.g. by choosing ( 2.11[ )(ii). A major advantage of (2.11) (ii) over (2.11 ) (iii) is 
that for the former it will be possible to derive almost linear finite element approximations 
that are unconditionally stable. The corresponding unconditionally stable schemes for the 
nonlinear shape function (2.11 ) (iii) , on the other hand, turn out to be more nonlinear. 
Conversely, if up > 0, then only g(— 1) = is needed in order to satisfy (2.16). The 



natural analogue for (2.11 ) (ii) in this situation is then 

g(s) = Ul + s), 



(2.17) 



and once again it is possible to derive almost linear finite element approximations that 
are unconditionally stable for this choice of g. 



Finally we note that the quartic potential (2.2) is often preferred in applications be- 
cause the discretized equations can then be solved with smooth solution methods, such 
as the Newton method. However, the quartic potential has the disadvantage that a priori 
it cannot be guaranteed that \tp\ < 1 at all times, and in practice it can in general be 
observed that discretizations of ip exceed the interval [—1,1]. Hence from a practical and 



from a numerical analysis point of view it is preferable to use the obstacle potential (2.3). 



Here we note that the discretized equations, which feature variational inequalities, can be 



efficiently solved with a variety of modern solution methods; see e.g. Barrett et al. (2004) 



Graser and Kornhuber 


( 


2007 


); 


Banas and Niirnberg 


(2008 


2009a 


); 


Blank et al. 


Hintermuller et al. 


(2011 


); 


Graser et al. 


(2012 


)• 





2011 



2.3 Anisotropies 



In this paper, we will only consider smooth and convex anisotropies, i.e. they satisfy 

j'(p)-q<j(q) G R d \ {0},q G M d , (2.18) 



which, on recalling (1.2), is equivalent to 



l(p)+i(p)-(q-p) <7(?) V P GM d \{0},gG 



(2.19) 



It is the aim of this paper to introduce unconditionally stable finite element approx- 
imations for the phase field models (2.7a-f) and (2.10a-d), (2.12a-c). Based on earlier 



work by the authors in the context of the parametric approximation of anisotropic geo- 
metric evolution equations Barrett et al. ( 2008afb ), the crucial idea here is to restrict the 
class of anisotropies under consideration. The special structure of the chosen anisotropies 
can then be exploited to develop discretizations that are stable without the need for 
regularization and without a restriction on the time step size. 

In particular, the class of anisotropies that we will consider in this paper is given by 

i 

L 



1 



(P)= X>€(P)F ), lt(p)-=\p-G t p]K WpeR d , re[l,oo), (2.20) 



where Gg G M dxd , for i = 1 — > L, are symmetric and positive definite matrices. This class 
of anisotropies has been previously considered by the authors in Barrett et al. (2008b 



2010b). We remark that anisotropies of the form (2.20) are always strictly convex norms. 



In particular, they satisfy (2.19). However, despite this seemingly restrictive choice, it 



is possible with (2.20) to model and approximate a wide variety of anisotropies that are 
relevant in materials science. For the sake of brevity, we refer to the exemplary Wulff 
shapes in the authors' previous papers |Barrett et al] ( |2008a| |b| |2010c|b|a| |2012a[ ). We 
remark that in the case r = 1 all of the numerical schemes introduced in Section |3j 
below, will feature no additional nonlinearities compared to the isotropic case (1.3). In 



particular, the finite element approximation in Section 3.1 for the obstacle potential (2.3 



will feature only linear equations and linear variational inequalities; see also Barrett et al. 



(2012c). Finally, we note that in the two-dimensional case (d = 2), the anisotropies (2.20) 
with the choice r = 1 adequately approximate most relevant anisotropies. However, in 
the three-dimensional setting (d = 3), it is often necessary to use r > 1 in (2.20) in order 



to model a chosen anisotropy. See|Barrett et al. (2008b) for more details 



In the following, we establish some crucial results for anisotropies of the form (2.20). 

7/(p)""' _i 



Note that for 7 satisfying (2.20) it holds that 

L 



A'(p) = 7 (p) y(p) , where Y(p) = ^ 



X L i(p) 



i t [p) VpeR d \{0}. (2.21) 



For later use we recall the elementary identity 

2y(y-z) = y 2 -z 2 + (y-z) 2 . 
Moreover, from now on we use the convention that 

liip) . * 



lip) 



if p = 



1 — >• L. 



(2.22) 



(2.23) 



Lemma. 2.2. Let 7 be of the form (2.20). Then it holds that 



7(p) < L'-c+d [Y^h 



r+l 



(2.24) 



.1=1 
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Moreover, 7 is convex and the anisotropic operator A satisfies 

A'(p) .(p-q)> 7 (p) h(p) - liq)} V P eR d \{0},qeR d , 

(2.25) 

hiq)]- 1 blip)? VpeR d ,qeR d \{0}, 

(2.26) 



A(p)<± 7 (<z)E 



lejp) 
.lip) 



r-l 



where in (2.26) we recall the convention (2.23) 



Proof. It follows from a Holder inequality that 

L 



r+l 



[j(p)y<l^ m 7 ^)r +i v P e 



which immediately yields the desired result (2.24). Next we prove (2.18). It follows from 



(2.21), a Cauchy-Schwarz and a Holder inequality that 

L r / \-ir-l I 



i'(p)-q = J2 



=1 



leip) 
lip) 



[llip)ViG,p).q<Y, 



= 1 



Itip) 
lip) 



-i r-l 



le(o) 



< 



7<(p) 
1 [lip) 



r-l 

r \ r / L 



5><W =7(9) V P GM d \{0},g 



Together with (1.2) this implies (2.19), i.e. 7 is convex. Multiplying (2.19) with 7 (p) 

Moreover, • 

^ b^p)Y 



yields the desired result (2.25). Moreover, we have from a Holder inequality that 

1 

L 



r+l 



[7(p)] r = I>(5)]^ 



< £b^)r £ 



[7<(p)] r ' 
L 

\lip)] r+1 < liq) $>^)] r+1 \lM)V VpeR d ,qe 



1 \ r+1 



This immediately yields the desired result (2.26), on recalling (2.5). □ 



Our aim now is to replace the highly nonlinear operator A'(p) : R — >■ R in (2.21) 



with an almost linear approximation (linear for r = 1) that still maintains the crucial 



monotonicity property (2.25). It turns out that a natural linearization is already given in 



(|2.21|). In particular, we let 

L 

7(?)E 



B r (q,p) := < 



7<(P) 
7(P) 



r-l 



~le(p} 
. liP) . 



r-l 



V J9 G 



(2.27) 



g = 
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where in the case p = we recall (2.23). For later use we note for q e M. d that 

B 1 (q,p) = B 1 (q,0) =: B^q) VpGM rf . (2.28) 

Clearly it holds that 

B r (p,p)p = A'(p) VpeR d \{0}, 
and it turns out that approximating A'(p) with B r (q,p)p maintains the monotonicity 



property (2.25). 



Lemma. 2.3. Let 7 be of the form (2.20). Then it holds that 



[B r (q,p)p].(p-q)>j{p) [7(p)-7(?)] Vp,gG 



(2.29) 



recalling (2.26), that 



Proof If p = then (2.29) trivially holds. Now let pel \ {0}. If q ^ it holds, on 



[B r (q,p)p] ■ (p ~ q) = l(q) 

L 

> 7(g) 

L r 

= 7(g) E 



Itip) 



-1 r-1 



[lM)\ (p-Q)-GtP 



1=1 



iip) 

nipY 



lip) 

nip) ibM)]' 1 iiip) - 1) 



iip) 



[nio)] 1 [nip)} 2 - iiq) iip) > iip) hip) - ii<i)} ■ 



If q = 0, on the other hand, then it follows from (2.24) that 



[B r (q,p)p] .(p-q) = [B r (0,p)p] .p = L*[7(p)] 1-r ^[7, 



r+1 



>[7(P)] 2 - 



t=i 



Corollary. 2.4. Lei 7 be of the form ( |2.20[ ). T/ien ZioZds t/iat 

[5 r (g, p) p].ip-q)> A(p) - A(q) Vp )9 eM d 



(2.30) 



Proof. The desired result follows immediately from Lemma 2.3 on noting the elemen- 



tary identity (2.22). □ 



3 Finite element approximations 

Let Q be a polyhedral domain and let {T h }h>o be a family of partitionings of Q into 
disjoint open simplices a with h a := diam(a) and h := max^-pi h a , so that Q = \J ae q-h~a. 
Associated with T h is the finite element space 

S h ■= { x e CiU) : x U is linear V cr G T' 1 } C F 1 ^). 
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Let J be the set of nodes of T h and {Pj}j£j the coordinates of these nodes. Let {Xj}j<=j 
be the standard basis functions for S h ; that is Xj £ S an d Xjipd = $ij f° r &U hj e ^- 
We introduce 7r h : C(fl) — >■ £> , the interpolation operator, such that (Tr h r])(pj] 
for all j G J. A discrete semi- inner product on C(fl) is then defined by 



v(Pj 



with the induced discrete semi-norm given by \r]\ h := [ ij],vi) h ]^ , for rj G C(O). We extend 
these definitions to functions that are piecewise continuous on T h in the usual way, i.e. 
by setting 

o-eT' 1 

where 

fc=0 



with {pj k }i =0 denoting the vertices of a, and where we define T)((pj k ) ) := lim r)(q), 
k = -> d. 

We introduce also 

# h := {x G 5 h : |x| < 1 in fi} C K := {r] G tf 1 ^) : \r)\ < 1 a.e. in fi} , 



7T 



{X G 5 h : X = on and S£ := { x e S h : X 

where in the definition of we allow for Ud G H*(dr>Q) fl C(dD&)- 



l u D on 9x)f2} , 



In addition to T h , let = t < t\ < . . . < tjy-i < t N = T he & partitioning of [0, T] 
into possibly variable time steps r n := t n — £ n _i, n = 1 — )■ iV. We set r := max n= i_»jv T n . 

In the following we will present stable finite element approximations for the phase 
field model (2.10a-d), (2.12a-c) for the obstacle potential (2.3) and for the case of a 



smooth potential such as (2.2), respectively. In order to obtain stable approximations, 



the three nonlinearities arising in (2.12a) from g(<p), from A'(V <p) and from ^/'(v 9 ) need 



to be discretized appropriately in time. Here the discretization of A'(V(p) induced by 



Corollary 2.4 is novel, and is one of the main contributions of this paper. The employed 
splitting of ^'(v 9 ) i 11 ^ implicit /explicit time discretizations according to a convex/concave 



splitting of on the other hand, is standard; see e.g. Elliott and Stuart (1993); Barrett 



et al. (1999). We employ the same idea to the splitting of g(ip), for which we now introduce 



some notation. A similar notation will be used in Section 3.2 for the splitting of ^'(v) m 
the case of a smooth potential 

Let G C 1 (1R) such that g(s) = g + (s) + g~(s). In our finite element schemes g + will 
play the role of the implicit part of the approximation of g, while g~ corresponds to the 
explicit part. We now define 



^(s ,si y 



g + (s 1 ) + g (s Q ) V s ,s x G 



(3.1) 
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as well as P ± (s) := J"* g ± (y) dy. Of particular interest will be splittings such that 

± u D (^)'(s) < V s < ^ . (3.2) 



If ud < 0, then (3.2) enforces P + (s) to be convex for s < -?=, while P (s) is concave over 
the same region. Possible splittings satisfying (3.2) for the shape functions in (2.11) are 
then given by 



(i) 




"00 


= 0, 




(«) 


= Q 


» 


l 

2 ' 




(ii) 






= 0, 


Q~ 


00 


= Q 


» 


= 1(1" 


- S 


iii) 


Q 




2 a ' 


Q~ 


00 


= Q 


» 


-1- = 


15 
16 



(3.3) 



(s 4 -2s 2 -|s + i; 



The fact that the splitting ( 3.3 ) (iii) satisfies ( |3.2 ) follows from the observation that in that 



case max < 2 g(s) 

•vs 



Q 



K 2 



-y/3< 



j| < |. Note that the above splittings were chosen 

such that the implicit part of the approximation of g is as simple as possible. If ub > 0, on 
the other hand, then swapping the roles of g in (3.3) will satisfy (3.2). However, as the 



implicit parts g + are then unnecessarily nonlinear in the cases (2.11 ) (ii) and (2.11) (iii), it 

(3.4) 



is more convenient, on recalling (2.17), to use the splittings 

1(1 + a), 



[11 Q' 



0, g (s) = g(s 
3 



g [s) 
g-(s) 



[ill) g'{ S ) = -^S 
which will then satisfy 

± Me^'OO < 



g{s) + f s 



i§( S 4 -2 S 2 + fs + l) 



V s > 



2 

V3 ■ 



(3.5) 



3.1 The obstacle potential 



We then consider the following fully practical finite element approximation for (2.10a-d) 



(2.12a-c) in the case of the obstacle potential (2.3). This approximation is an adaptation 

can 



of the scheme from Elliott and Gardiner (1996) which, with the help of Corollary 2.4 
be shown to be stable. Let $ u G K h be an approximation of (po G K, e.g. $° = w^fa for 
ipo G C(Q). Similarly, if •& > let W° G be an approximation of Uq. Then, for n > 1, 
find ($ n , W n ) eK h x S h D such that 



'W n -W n ~ l \" / , $ r 



TV] 



71-1 



TV) 



..V) +(7r ft [6($^)]V^,Vx) = 
V x e Sft, (3.6a) 



/i(V $ n_1 ) — , x - $ n ) + ^ (S r (V $ n_1 , V $ n ) V $ n , V [x - $ n ]) 

> (c*^g(<5> n - 1 ,<5> n )W n + e- 1 <5> n - 1 ,x-$ n ) h W X ^ K h . 



(3.6b) 



The main differences between (3.6a,b) for g + = 0, so that £)($" 1 , $ n ) = g(<$> n 1 ), and the 



basic scheme in Elliott and Gardiner (1996, Eqs. (3.1), (3.2)) are our novel approximation 
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of A'(V <p) in (3.6b) and the fact that we evaluate the discrete temperature on the new 



time level W n in (3.6b). The latter implies that the system (3.6a,b) is coupled, and this 

below. We stress that 



is needed in order to derive a stability bound, see Theorem 



3.6 



there is no stability result for the scheme Elliott and Gardiner (1996, Eqs. (3.1), (3.2)). 



In addition, we allow for the splitting g = g + + g , so that unconditional stability can 
still be shown for nonlinear functions g. 



Let 



x) 



\W-u D \l + 



Xa 1 
a c^j 



[feKV^ + e- 1 (¥(<&), l) h ] 



and define 



J*(W, $) = £~(W, $) - Xu D (P($), 1)' 



for all W, G S h , as the natural discrete analogue of the energy appearing in (2.13). We 
can then show that the solutions to (3.6a,b) satisfy a discrete analogue of (2.13). 



We begin with considerations for the almost linear scheme (3.6a b) with g 
-- 1. 



1, let g 



■ and 
and let ud G M. 



Lemma. 3.1. Let 7 be of the form ( |2.20[ ) with r 
Then there exists a solution W n ) G K h x to (3.6a,b) and $ n ; W n are unique 
up to additive constants. If p + (|o(<I> ri-1 )|, l) h + |($ n ~ , 1)| > 0, then <£> n zs unique. If 
d > or c^Affi 7^ <9f2 t/ien zs unique if $ n unique. If r) = and 8nQ = dQ, and 
if $ n is unique, then W n is unique if there exists a j G J such that \$ n (pj)\ < 1 and 
g^-HPi)) ^ 0. 



Proof. The proof follows the ideas in Barrett et al. (1999), see also Blowey and Elliott 
fll992|). At first we assume that d N Vt ^ 5fi or that ■& > 0, so that £ h : S' 1 such that 



(vr^^ 1 )] V [Q h v%Vrj) + - {g h v\rj) h = (v h ,r]) h V 77 G S* , 

is clearly well-defined, on recalling that 

b(s) > min{/C+,/C_} > VsG [-1,1]. 
Moreover, it follows from ( |3.6a[ ) and (3.7) that 



v h G S h (3.7) 



w n -u D = g h 



(W 71 ' 1 -u D ) - Xn h 



\n-l 



71-1 



(3. 



Substituting g§ into ( |3.6b[ ), and noting Q with w' 1 = vr /l [^(<l>"- 1 , $ n ) ( X - $")] and 
77 = £ ft n h [g{®^, $ n ) ($ n - $ n ~ 1 )] yields that 

Xc^j a 



> /l [6($ n - 1 )] V \Q h 7T h [g(^ n -\ $ n ) ($ n - $™~ 1 )]], V 7r fe [^($ n - 1 , $ n ) ( X - $ r 



1? 



(^ 7 r' l [^($ ri - 1 ,$ n ) ($ n - $ n - 1 )],^7r /l [^(<l> n - 1 ,<I) n ) (x - $ n )]) /l } 



+ (/x(V $ n , x - $ n f + e {B r {V $ n -\ V $ n ) V V [x - $ n ]) 



CtTr, 
eh 



>(Ax-$ 



n\h 



\/ x eK h 



(3.9a) 
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where 



ep 



/i(V$ 



n-l 



n-l 



+ — h 



n-l 



«o]) 

(3.9b) 

is piecewise continuous on T h . As we consider the case g + = 0, from now on we use the 
fact that ^($ n - 1 ,$ n ) = We recall from fl2.27| ) an d (j2 28ft that 5i(g) £ M dxd is 

symmetric and positive definite for all q £ M d , and hence (3.9a) are the Euler-Lagrange 
equations for the convex minimization problem 



mm 



A &j, a 
2ar n 



x) 



2r„ 



(7r /l [6($ n - 1 )], |V [G h ir h { g (<f> 



+ 



2ar„. 



(m(v$ 



2\/i 



+ -( J B 1 (V$"- 1 )Vx,Vx)-(/\x) / 



Therefore there exists a $ n £ If' 1 solving (3.9a) that is unique if p > or 7r /l [^($ n ~ 1 )] 7^ 
£ S h , and is unique up to an additive constant otherwise. In the latter case, if 
1) 7^ 0, then it immediately follows from (3.6b) that $ n is unique. If $ n is unique, 
then the existence of a unique W n £ Sp, such that ($™, W™) solve (3.6a,b), follows from 

Q. 

For the remainder of the proof we assume that c^fi = <9f2 and that 1? = 0. Then it 
follows immediately on choosing x = 1 in (3.6a) that (£>(<3> n_1 ), = (^(J?™ -1 ), $ n ~ 1 ) /l . 
Taking this into account, we define Q h : — >■ S* such that 

(tt 71 ^"- 1 )] V[Q h v h \,Vrj) = (v h ,T]) h \/rieS\ v h £ S h , 

where S h := {x £ S'' 1 : (x, 1) = 0}, and observe that (3.6a) then implies that 

A 



g h 7i h [gi^- 1 ) ($ n - $ n " 1 )] + c 



where £ n £ R is a Lagrange multiplier. It follows that (3.9a) holds with r) 
replaced by Q h , and with K h replaced by K 



(3.10) 

0, with g h 



{ X eK h : ( Q (<S> n - 1 ), X -$ n - 1 ) h = 0}- 
As before we can interpret this variational inequality as the Euler-Lagrange equations of 
a convex minimization problem, which yields the existence of a solution $ n £ K h that 
is unique unless p = 0, 7T h [g(^ n ^ 1 )] = and ($ n_1 ,l) = 0. Therefore, on noting (3.10), 
we have existence of a solution ($ n , W n ) £ K h x S h to (3.6a b). If $ n is unique, and 
if |$ n (pj)| < 1 and ^($ n_1 (pj)) 7^ for some j £ J then (3.6b) holds with equality for 
X — Xji which uniquely determines £ n and hence yields the uniqueness of W n . □ 



It turns out that most of the technical assumptions in Lemma 3.1 are trivially satisfied 
for the shape function choices in (2.11). In particular, we obtain the following result. 



Corollary. 3.2. Let 7 be of the form (2.20) with r = 1, let g be given by one of 
the choices in (2.11) or by (2.17), let g + = and let up £ M. Then there exists a 
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solution ($ n , W n ) G K x S D to (3.6a,b) and $ n , W n are unique up to additive constants. 
Moreover, <3> n is unique unless g is of the form (2.11)(iii) ; and p = 0, ( | <^> ri ~ 1 1 , \) h = \Q\ 
and =0. 



Proof. The desired results follow immediately from Lemma 3.1 □ 



Remark. 3.3. Let the assumptions of Lemma |3.1 hold and let d^Q ^ dQ. Then it is 



easy to prove that if $ n 1 = 1 and $ (W n 1 — ud) = 0, and if 



oil) u D < — e 1 



(3.11) 



then the unique solution to (3.6a b) is given by $ n = 1 and W n = ud- If the phase field 
parameter e does not satisfy (3.11), then $ n = 1 and W n = ud is no longer the solution 



to (3.6a,b). In practice it is observed that if e does not satisfy (3.11), then the solution $ n 



exhibits a boundary layer close to dVt where $ n < 1. This artificial boundary layer is an 



undesired effect of the phase field approximation for the sharp interface problem (1.4a-e 



In fact, and not surprisingly, (3.11) is precisely the condition on g(l) in (2.16). This 



motivates the use of shape functions with g(l) = 0, such as (2.11 )(ii) and (2.11)(iii), in 
practice. An obvious advantage over e.g. (2.11 )(i) then is to be able to use larger values 



of e, which in itself means that less fine discretization parameters may be employed. 



For completeness we note that if, and only if, the condition 



a 



l)u D < —e 



-i 



holds, then $ n 
(W n ~ l 



•1, W n = ud is the unique solution to (3.6a,b) for $ 



Ud) = 0. Satisfying both (3.11) and (3.12) is equivalent to satisfying (2.16) 



\n— 1 



(3.12) 
-1 and 



Remark. 3.4. Let 7 be of the form (2.20) with r > 1, and let the remaining assumptions 
of Lemma 3.1 hold. Then the highly nonlinear system (3.6a b) for ($ n , W n ) is no longer 



continuously dependent on the variable $ n , recall (2.27). Due to this fact it is not possible 



to show existence of solutions to (3.6ab) with the help of Brouwer's fixed point theorem. 



However, in practice we have no difficulties in finding solutions to the nonlinear system 
(3.6a b), and the employed iterative solvers always converge; see Section 4.2, We recall 



that the same situation occurred in Barrett et al. (2008b), see Remark 3.3 there, where 



discretizations for anisotropic geometric evolution equations for anisotropic energies of 



the form (2.20) were considered for the very first time. 



We now extend the existence result from Lemma |3.1| to the case of a general splitting 
g = g + + g~ . On recalling from (2.16) and from Remark 3.3 that nontrivial choices of g, 
i.e. alternatives to (2.11 )(i), are only of interest when d^Q 7^ dQ, we consider the case 
g + 7^ only in the presence of Dirichlet boundary conditions on W n . 



Theorem. 3.5. Let 7 be of the form (2.20) with r = 1 and let ud £ R- In addition let 
p+ |($ n -\l)| > or 

(|^- 1 , X )|,l) ft >0 V X eK h . (3.13) 
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Moreover we assume that either g + = 0, or d > 0, or d^fl ^ dQ. Then there exists a 



solution (<f> n ,W n ) eK h xS h D to (3.6a,b) 



Proof. The desired result for the case g + = has been shown in Lemma 3T We 
now consider the case g + ^ 0, so that either d > or d^fl ^ dQ. Then we can apply 
Brouwer's fixed point theorem to prove existence of a solution $ n as follows. Let the 



map T : K h -> K h be defined such that $ new = T($ old ) is the solution of ( |3.9a[ b) with 
^($ n_1 ,$ n ) replaced by ^"($ n_1 , $ old ), and with all other occurrences of $ n replaced by 



$ ncw . It follows from the proof of Lemma 3A_ and our assumptions that there exists a 
unique $ now G K h , and the continuity of the map $ old h-> $ ncw = T($ old ) together with 
the fact that K h is compact and convex then yields the existence of a solution $ n G K 



to (3.6a b). The existence of a solution W n G then follows from (3.8). □ 



The following stability theorem is the main result of this paper. 



Theorem. 3.6. Let 7 be of the form (2.20) and let ud G KL Then it holds for a solution 
($» w n ) eK h x S h D to (|&6a|b) tfmi 



£ T A (r, $ n ) - w D A (g^* -1 , $ n ), $ n - $ n ~ 1 )' 1 + r n (7r /l [6($"- 1 )] V W", V W n ) 



Xp € 

a cq, 



[MVf- 1 )] 



1 $ n - $ 



n-l 



n— 1 /Rn— !"> 



(3.14) 



In particular, if the splitting g = g + + g satisfies 

±u D (g ± )'(s)<0 VsG[-l,l] 

then it holds that 

Xp e 



(3.15) 



^(W n , $ n ) + r n (7r /l [6($ n - 1 )] V W n , V W n ) + r n 



Kv^- 1 )]^ 



Proof. Choosing \ — W n — ud in (3.6a) and \ = & n i n (3.6b) yields that 



7 ' 

n-l 



(3.16) 



+ r n (7T h [6($ n - 1 )] V iy n , Vr) = 0, (3.17a) 
e ^ r- 1 (/i(V $ n_1 ) $ n - $ n_1 , $ n " 1 -$ n ) h + e (B r (V $ n ~\ V $") V $ n , V [$ n_1 - $ n ]) 

> (c* ^ £?($ n -\ $ n ) + e" 1 $ n_1 , $ n_1 - $ n )\ (3.17b) 



a 



It follows from (3.17a,b), on recalling (2.22) and (2.30), that 

■& a 



|e|7(V$ 



n\ 1 2 1 _— 1 I ^f>n 1 2 



$"|2 + 



2 , _ _ P 



2 Aa 



c^\W n -u D \i +r n e 



a 



1 $ n - $ 



n— 1 



h 



A a 



<i£| T (V$ 



n-lM2 1 _-l |^n-l|2 



7? a 
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This yields the desired result (3.14) on adding the constant \ e 1 f n 1 dx on both sides, 



and then multiplying the inequality with — . In addition, it follows from $ n , $ n £ X 



and (3.15) that 



< w D (p-($ n ) - p-($"- 1 ) + P+($ n ) - P+($"- 1 ), if 
= u D (P($ n ) -P($ n - 1 ),l) ft . (3.18) 



The desired result (3.16) now follows on applying (3.18) to (3.14). □ 



3.2 Smooth potentials 



The unconditionally stable approximation (3.6a,b) for the obstacle potential (2.3) can be 
easily adapted to the case of a smooth potential such as (2.2). To this end, let := 
for an arbitrary smooth potential and let = <p + + 0~, with <f> being the derivatives of 
the convex/concave parts of i.e. 



± (0 ± ) / (s) > V s £ 



(3.19a) 



We will make the mild assumption that there exist constants i/)o,ipx,8 > such that 

V + (s) > *pi \s\ 1+s -ip VsGl. (3.19b) 



For the quartic potential (2.2) the natural choices are 

4> + (s) = s 3 and 



s = — s . 



(3.20) 



so that (3.19a,b) are clearly satisfied. 



As before, given $° £ K h and, if i? > 0, W° £ for n > 1, find ($ n , W n ) G5 ft x^ 
such that 



n-l \ 71 



n-l \ ft 



- (7T /i [6(<i> n ~ 1 )] vr,vx) = o 

V X e S \ (3.21a) 



e- I u(V$ 



1 $n _ fon-l 

n— 1\ 



, X ) + e (S r (V $ n ~\ V $") V $™, V X ) + e~ l (</>+($"), x) h 



ai-l\ 



where in order to avoid degeneracies we have defined 

[6(1) s >!> 
6(s) = < b(s) \s\<l, 

Ife(-l) s<-l 



(3.21b) 
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and where for technical reasons we have introduced 

o + (m) s > m . 

Qm(so,Si) ■= Q~(s ) + £>+0i) , where £+(s) := { g + (s) \s\ < m, (3.22) 

Q + (—m) s < — m . 



for some fixed parameter m > 2. We note that these modifications of (2.8) and (3.1) are 
such that 



and 



max \Qm{so,s 



b(s) > min{/C+,/C_} > 

C(m,s ) 



= max \g(s , s ) 

\s\<m 



V seR, 
Vs Gl, 



(3.23a) 
(3.23b) 



Theorem. 3.7. Let'j be of the form ( |2.20| withr = 1 and letu D eR. If q + = andif4> + 
is strictly monotonically increasing, then there exists a unique solution (<3> n , W n ) G S h xS f l) 



to (3.21a,b). If q + 7^ 0, and if either $ > or d^jVt ^ dVt, then there exists a solution 
G S* x S£ to (|3.21a[b) if ^+ satisfies the assumption fl3.19b[). 



Proof. The existence and uniqueness proof for the case g + = 0, which is a simple 
modification of the proof of Lemma 3A_, is left to the reader. Note that this proof makes 
use of the strict monotonicity of + . 

In order to proof existence for the case g + ^ 0, we apply Brouwer's fixed point theorem. 
It is this part of the proof that requires the cut-off of ^ defined in (3.22), as well as the 
mild assumption (3.19b). The application of Brouwer's fixed point theorem is similar to 
the proof of Theorem 3.5 Setting up the map 

$°id ^ $new = T ( $ oid) ana i g 0US i y to the 
proof there, we immediately see that the map T is well-defined and continuous, where 
we recall that our assumptions yield that $ > or ^ dQ. It remains to show that 
T : Y h — > Y h for a bounded subset Y h C S h . To this end, on recalling (3.9a b), we note 
that $ ncw G S h satisfies 

A cq, a 



n h [b($ n - 1 )] V \Q h 7r h [Q m ($ n -\ $ old ) ($ new - $ n_1 )]], V \Q h TT h [g m (<S> n -\ $ old ) x])) 



n — X ,jjold\ /^ncw 



^ 1 )],0 ft 7r ft [£ m ($^\$ QU )x]) 



+ 



ep 



ar n 

-h Aft 



(Ai(V$ 



n-l^ 



X) h + e (^(V & n ~ L ) V $ new , V x) + e- 1 (0 + ($ new ), x) 



(g,xY 



V x e S* 



where g h :- 



Ai(V$ 



n-X 



rn-X 



ud\), and where Q h is defined by (3.7) with b replaced by 6. These are the Euler-Lagrange 
equations for the convex minimization problem 



min [J( x )-(g h ,x) h ] 



(3.24a) 
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where 



J(x) 



Ac$a 
2ar n 

+ 



t) 



2r, 



(-K h [b($ n -%\V[g h -K h [Qm(® 



"- 1 $ old ) _ $ n ~ l 



;jji 2 ) 



+ - (B 1 (V V X , V X ) + £ _1 (^ + (x), 1 



2 a r, 

h V x e S 71 



2\/i 



It follows from (3.24a,b) and (3.23a b) that 



1 (^+($ r 



(/.r")" < J(0) < C($ 



Applying the elementary inequality y z < ^ \y\ p + i for p,q & (1, oo) with ^ 
to the second term in (3.25) yields that 



(3.24b) 
(3.25) 

(3.26) 



where and 5 are as in (3.19b). Now combining (3.25) and (3.26), on recalling the mild 
assumption ( 3.19b ), yields that (|$ new | 1+5 , l) h < C for some constant C > independent 
of $ old , i.e. 

(|T( X )| 1+ * l) h <C7 \/ X eS h . (3.27) 

Hence T : Y h — > Y h for a bounded subset Y h C S h , and so Brouwer's fixed point theorem 
yields the existenc e of a solution $™ to ( 3.21a b). The existence of a solution W n G Sy^ 
then follows from (3.8) with Q h replaced by Q h and with ^"replaced by m . □ 



We stress that the cut-off introduced in (3.22 ) is for technical reasons only. If a solution 
($ n , W n ) G S h x S h D to fl3 21a| b) is such that |$ n | < m, then clearly ($ n , W n ) G S h x S h D 
also solves (3.21a,b) with ~g m replaced by q. In practice, this is always the case for m 
chosen sufficiently large. Hence for practical implementations, only ( 3.21a[ b) with £ m 
replaced by g needs to be considered. 



Corollary. 3.8. Let 7 be of the form (2.20) with r — 1 and let u D e R. Let \I> be given 
by ( 2.2[) and Ze£ (3.20) ZioZd. Lei £ and zts splitting be defined by one of the choices in (3.3) 
or (3.4). Then there exists a solution ($ n , W n ) £ S h x S^, to (3.21a,b) unless in cases 
(3.3) (iii) and (3.4)(iii) it holds that ■& = and d^Cl = OVl. Moreover, ($ n ,W n ) is unique 
for the choices ( 3.3 ) (i) , ( 3.3 ) (ii) and ( 3.4 ) (ii) . 



Proof. The desired results follow immediately from Theorem 3.7 on noting that (3.20) 
satisfies the assumptions on <p + and \E' + stated there. □ 



(3.21a,b) when r 
z f = 

$ n = 1 and W r 



Remark. 3.9. Similarly to Remark 3.3, the following observation holds for the scheme 
1; Q + — 0? ud G R and <p + is strictly monotonically increasing. Then, 
and "d {W n ~ l — Up) = 0, then the unique solution to (3.21a,b) is given by 
up if and only if 



udq{1) = 0. 



(3.28) 
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For nonzero up this is precisely the necessary condition for G(s) in (2.15) to have a local 
minimum at s = 1. In practice, if the condition (3.28) is violated, then for certain values 



ofuo ande artificial boundary layers develop. This undesired effect for the choice ( 2.11[ )(i 
once again motivates the use of the alternatives (2.11 )(ii) and (2.11) (iii) in practice. 



The following stability result is the natural analogue of Theorem 3.6 for the case of a 
smooth potential 



Theorem. 3.10. Let'y be of th e form ( 2.20) and letuo G H£. Then it holds that a solution 
($ n ,W n ) e S h x to (3.21a b) satisfies (3.14) with b replaced by b, and with g replaced 



by g m . In particular, if the splitting g = g + + g satisfies (|3.2|) ; and if 

2 

V3 



n— 1 



<7I and 



or if it satisfies (3.5), and if 



n-l 



> 



(3.29a) 



and $" > 



V3 



(3.29b) 



then the solution ($ n , W n ) satisfies the stability bound (3.16) with b replaced by b. 



Proof. The proof of the stability bounds, which is a simple modification of the proof 
is left to the reader, 
and g = g + + g~ 



of Theorem |3.6 

= </> + + 



Note that the proof makes use of the splittings 



recall (3.18). □ 



Corollary. 3.11. Let 7 be of the form (|2.20|) and let u D e 

/ • 

1 



of g and its splittings in (3.3) 
($ n , W n ) eS h x S h D to J3.21a 



Then for the choices 
( 3.3 ) (ii) and (3.4 ) (ii) it hol ds tha t the unique solution 
b) satisfies the stability bound (3.16) with b replaced by b. 



For the choice (2.11 ) (iii) ? with the splittings ( 3.3 ) (iii) or (3.4)(iii), it holds that a solution 
($» W n ) G S h x S h D to ( |3 21a| b) satisfies the same stability bound if ( |3.29a[ ) or ( |3.29b[ ) 
hold, respectively. 



Proof. The desired results for the splittings (3.3)(i), (3.3 ) (ii) and (3.4 ) (ii) , on recalling 



Corollary 3.8, follow from the fact that these splittings satisfy the inequalities in (3.2) for 
all seR. The results for (2.11 ) (iii) follow immediately from Theorem 3.10 □ 



In practice, in general, the values of $ n are either within the interval [—1, 1], or very 



close to it. In our experience, for the scheme (3.21a,b) with (2.11 ) (iii) and with the 



splitt ings ( 3.3) (iii) and ( 3.4 ) (iii) for ud < and «d > 0, respectively, in practice (3.29a) 
and (3.29b) always hold. Here we note that -7= ~ 1.15. 



4 Solution of the algebraic systems of equations 



The system of nonlinear equations for ($ n , W n ) arising at each time level from the ap- 



see e.g. 



proximation (3.21ab) can be solved with a Newton method or with a nonlinear multigrid 
method, 



Kim et al. (2004) 
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For the remainder of this section we discuss the solution of the systems of algebraic 
equations for ($ n , W n ) arising at each time level from the approximation (3.6a b). Adopt- 
ing the obvious notation, the system (3.6a b) can be rewritten as: Find ($ n , W n ) G 
[-1, 1] J x R J , J := #J, such that 



A M e ($ n ) $ n + (t? M + r n A) W r 



P 



e-r- 1 (V-$ 
a 



n\T 



+ e(V 



- /($ n ) 
$ n ) T B r ($ r 



(4.1a) 



$ n - - (V - $ 



> (V- $ n ) T # 



V 7 G [ 



(4.1b) 

where M, M M , M e (r]), A and B r {rj), for 77 G S'' 8 , are symmetric J ^ J matrices. In the case 
of pure Neumann boundary conditions, (1.5)(ii), their entries are given by My := (xi, Xj) h i 

V) Xu Xj) h , 



[M, 



n—l 



while the right hand sides in this case are defined as /($ n ) := A M £ ,($ n ) $ 



n— 1 



and g 



n-l 



( 1.5[ )(iii) these entries need to be appropriately manipulated. 



Of course, for the cases ( 1.5 ) (i) and 



n-l 



Clearly, the algebraic system (4. lab) can be written as a (symmetric) nonsmooth 



saddle point problem of the form: Find (U, W) G [—1, 1]^ x ¥L J , 

Mg(U) U + AW = fg(U) 

(V - U) T C r {U) U - (V - U) T M e (U) W>(V-U) T g V V G [-1, 1] 



J 



(4.2a) 
(4.2b) 



where we prefer to write the unknowns as (U, W) in place of ($ n , W n ), in order to highlight 
the connection to discretizations of Cahn-Hilliard equations, where the former notation 
is standard. On recalling (2.28) and (3.1), we note that (4.2a,b) in the case r = 1 and 
q + = collapses to 



MU+AW = f 



(V -Uf'CU - (V -U) T MW > (V~U) T g VVg [-1,1] 



J 



(4.3a) 
(4.3b) 



where C := Ci(0), M. := A^ e (0) and / := f e (0). Nonsmooth saddle point problems of 
the form (4.3a,b) are well-known from the numerical approximation of (isotropic) Cahn- 
Hilliard equations. Various different solution methods for the system (4.3a b) are discussed 
inlBarrett et al. I (120041) ;|Graser and Kornhuberl (120071) ;|Banas and Niirnbergl (120081 |2009a|); 



Blank et al. (2011); Hintermuller et al. (2011); 



Griiser et al. (2012). In the case r = 1 we 



use the solution method from Banas and Niirnberg (2008) in order to solve (4.3a b). In 
the remainder of this section we consider the case r > 1. We now state possible solution 
methods for the nonlinear nonsmooth saddle point problem (4.2a,b). 



4.1 Nonlinear Uzawa-multigrid iteration 



In what follows, we will extend the Uzawa-multigrid iteration from Bahas and Niirnberg 



(2008), which is based on the ideas in Graser and Kornhuber (2007), to the highly non- 



linear saddle point problem (4.2a b). The method from Banas and Niirnberg (2008) can 
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be interpreted as a primal active set method, where the approximation of the active set 



Blank et al. (2011) 



is driven by the current iterate W k in (4.3b), rather than via a dual parameter as in e.g. 



Given an initial iterate (U Q , W ) G [-1, 1] 3 x R J , for k > let U k+ i G [-1, 1] J be the 
solution of 



(V - U k+l fC r {U k ) U k+h >(V- U k+h ) T (g + M e (U k ) W k ) V V G [-1, if . (4.4a) 



Then we define the active sets as 

Jt ■= {j G J : [U^i]j = ±1} and let J k = J+ U J k . 
Now we seek the solution (U k+ i, W k+ i) G M. J x M. J to the linear system 

C r (J k ,U k ) -M e (J k ,U k )\ ( U k+1 \ = fg(4,J k 
M e {U k ) A \W k+1 { f e (U k ) 



(4.4b) 



(4.4c) 




i G J k , 



[C r (U k )] l3 ieJ\J k 



[M g (J k ,U k )]i 



±1 zeJ±, 

gi i G J\J k 



i G J k , 



Now we continue the iteration (4.4a-c), until convergence is obtained, i.e. until 



J, 



± 

k+l 



Jt 



and max { max | [C/jfe+i]j 



[U k }j \ ,max|[WAH 



< tol 



(4.5) 



where to/ 
then for k 



10 8 is a given fixed tolerance. If a good initial guess Wo is not available, 
it can be beneficial to set Ui 



Uo, rather than to employ (4.4a). Observe 
that since the iterates U k+ i are only needed to define the active sets J k in (4.4b), an 



iterative procedure to find the solution of (4.4a) can be stopped as soon as the active sets 



J k have been found. In practice we stop the iteration as soon as two successive iterates 



for (4.4a) have the same active sets, which is usually the case after a few projected block 



Gauss-Seidel iterations. Alternatively, a monotone multigrid method could be employed 



to solve (4.4a), see Kornhuber (1994). The linear saddle point problems (4.4c) can be 



solved with a multigrid method using block Gauss-Seidel smoothers or, alternatively, 



with a direct solution method such as UMFPACK QDavis| (|2004[)) or LDL QDavis| Q2005fl ), 
together with the sparse matrix ordering package AMD (Amestoy et al. ( 2004[)). H ere for 
the multigrid solver and the LDL factorization package, the linear system (4.4c) needs 



to be equivalently reformulated with a symmetric block matrix, which is easily possible. 
Finally, we observe that in the case r = 1 and q + = 0, the first stopping criterion in (4.5) 



immediately implies the second criterion in (4.5), as then the linear system (4.4c) does 
not depend on the iterates U k . 
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Remark. 4.1. In practice, in our computations, the iteration (4.4a -c) did not converge 
for values of r > 3, while it usually converged for smaller values of r. In particular, 



it always converged in the case r = 1 for the nonlinear approximation (3.6a b) with the 



splitting (3.3)(iii). However, as we are interested in performing simulations for much 



larger values of r, e.g. r = 9 for ani 9; below, we also consider a more robust solution 
method in the next subsection. 



4.2 Lagged fixed point iteration 

In this subsection we consider a lagged fixed point iteration, where at each iteration a 



subproblem of the form (4.3a b) needs to be solved 



Let k = 0. Given an initial iterate (U , W ) G [-1, 1] J x R J , we let (U k+ i,W k+ i) G 
-1, l]" 7 G xR J be the solution of 



M e {U k )U k+l + AW k+k = f e {U k ) 
(V - U k+h ) T C r (U k ) U k+h -(V- U k+h ) T M g (U k ) W k+h > (V — U k+h ) T g 



(4.6a) 



V V G [~l,l} J ■ (4.6b) 

On obtaining (U k+ i,W k+ ±), we set 

(U k+1 , W k+1 ) = (1 - //) (U k , W k ) + /i {U k+h , W k+ i) , (4.6c) 



where ji G (0, 1] is a fixed relaxation parameter. The iteration (4.6a -c) is repeated until 



max \ max|[C7 fe+ i]j- - [U k ]j \ ,max\[W k +i]j - [W k ]j\ > < tol . 



In practice, the iteration (4.6a -c) always converged, provided \x was chosen sufficiently 
small. 



5 Numerical experiments 



In this section we report on numerical experiments for the proposed finite element approx- 



imations. Apart from a single computation for the approximation ( 3.21a b) in the case of 



the quartic potential (2.2), where we employ the splitting (3.20), we will present results 



for the approximation (3.6a b) for the obstacle potential (2.3) only. Our preference for 



the scheme (3.6ab) over the alternative approximation ( 3.21a b) stems from the fact that 



in the former the phase field approximation $ n is guaranteed to stay inside the interval 
[—1,1], while the latter scheme in general admits values |$ n | > 1, which in practice is 
observed if e.g. a well developed interface is present. Moreover, the bulk regions for the 
approximation (3.6a,b) are easily identified through $ n = ±1, whereas for the scheme 
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(3.21a,b) this is less straightforward. For the implementation of the approximations we 



have used the adaptive finite element toolbox ALBERTA, see Schmidt and Siebert (2005) 



For the approximation (3.6ab) we employ the adaptive mesh strategy introduced in Bar- 
rett et al. (2004) and Banas and Niirnberg (2008), respectively, for d = 2 and d = 3. This 
results in a fine mesh of uniform mesh size hf inside the interfacial region |$ n_1 | < 1 
and a coarse mesh of uniform mesh size h c further away from it. Here hf = and 

h c = ^ are given by two integer numbers Nf > N c , where we assume from now on that 
Q = (-H, H) d . In all of the experiments below we have H = 
otherwise stated. 



with ( 1.5 ) (i) , unless 



Throughout this section the initial data (p G C(fi) is either chosen constant, ip = 1, 
or is chosen with a well developed interface of width e it, in which ip Q varies smoothly and 
such that r = {x G Q : ipo(x) = 0}. Details of such initial data can be found in e.g. 
Barrett et al. ( |2004 ); Banas and Niirnberg (2009b, 2008). In general the initial interface 
r is a circle/sphere of radius Ro G (0,H) around the origin. We use Rq = 0.1 unless 
otherwise stated. If $ > 0, we set 







u {z) 



1 - e R( > 
u D 



\z\ < Ro , 
— (1 - e«o-M) R <\z\<H. 

\z\ > H. 



We always fix $° = n h ip and, if i? > 0, W c 



7T h U Q . 



Unless otherwise stated we always let e^ 1 = 16 it and Nf = 128, N c = 16. In addition, 
we employ uniform time steps r n = r, n = 1 — > N. As an indication for the computational 
effort that is involved in producing the simulations presented in this section, we state for 
each simulation an exemplary CPU time for a single-thread run on an Intel i7-860 (2.8 
GHz) processor. 

For the anisotropics in our numerical results we always choose among 



ANI 



(*). 



y(p) = J2 [5 2 \p\ 2 +p)(l-5 2 )Y , with 5>0 

3=1 



ani 2 : 7 as on the bottom of Figure 3 in Barrett et al. (2008b) 



ANI3: 7 as on the right of Figure 2 in Barrett et al. (2010b), 



ANI4: 7 as in Figure 3 in Barrett et al. (2012b) , 

ANlg: 7 as on the right of Figure 3 in Barrett et al. (2010b) 



r = 


1 


,L 


= d], 


r = 


1 


,L 


= 2], 


r = 


1 


,L 


= 3], 


r = 


1 


,L 


= 4], 


r = 


9 


,L 


= 3]. 



We remark that ani^ is a regularized Zx-norm, so that its Wulff shape for 5 small is given 
by a smoothed square (in 2d) or a smoothed cube (in 3d) with nearly flat sides/facets. 
Anisotropies with such flat sides or facets are called crystalline. Also the choices ANI,, 
i = 2 — > 4, represent nearly crystalline anisotropies. Here the Wulff shapes are given by 
a smoothed cylinder, a smoothed hexagon and a smoothed hexagonal prism, respectively. 
The Wulff shape for the cubic anisotropy ANlg is given by a smoothed octahedron. Finally, 
we denote by ANl£ the anisotropies ANI&, k — 3 — > 4, rotated by i n the X\ — x 2 -plane. 



26 




Figure 1: (e^ 1 = 16 tt, ani!°- 01) , (|2.1lfc(i), a = 1, p = 10~ 3 , u D = -65, Q 



Creation of a boundary layer for the scheme (3.6a,b). Snapshots of the solution at times 
0, 4 x 10" 5 , 5 x 10" 5 , 7 x 10" " 



1 

2' 2' 



t 



10" 



Finally, unless otherwise stated, we choose A 



fC± = 1 and (3 = 7, where we 



recall (2.6). 



5.1 Mullins— Sekerka in two space dimensions 



In this subsection we always choose 1? = 0. We begin with an investigation into the choice 
of q. At first we choose (2.11 )(i). In order to visualize the possible onset of a boundary 



layer as explained in Remark 3.3 for the obstacle potential (2.3), we present a computation 
for (3.6a,b) with the initial data $° = <^ = 1. For this experiment we use p 

= -± 16tt 



1, the critical value for ud in (3.11) is — ^-e 



lO" 3 . On 
—64. In our 



setting a 

numerical computations this lower bound appears to be sharp. In particular, we observe 
that U n = 1 is a steady state whenever «d > —64, but a boundary layer forms already for 
e.g. Ud = —64 — 10~ 8 . The same behaviour has been observed by the authors in Barrett 
et al. (2012c) for the choice p = 0. As an example for the case p = 10 -3 considered here, 



we present a run for ud = — 65 in Figure [TJ where we can clearly see how the boundary 
layer develops. Note that this phenomenon is completely independent from the choice of 
anisotropy 7. The discretization parameters for this experiment were Nf = N c = 128 and 
t = 10~ 5 . Note that in the presence of p > we observe a convex shape in Figure [TJ in 
contrast to the corresponding evolution in Barrett et al. (2012c, Fig. 4), where p = 0. 



For the approximation (3.21a b), i.e. in the case of the smooth quartic potential (2.2) 



we observe that the criterion ( 3.28[ ) is of course not sharp, in the sense that even for 
values > ud > —41 no boundary layer forms in practice, even though (3.28) is then 



violated for the shape function (2.11 ) (i 
values less than 1 



What happens in practice is that $ n attains 
without forming an interface, i.e. min xg n <& n {x) > 0. However, for 



the value ud = —42, with the remaining parameters fixed as in Figure [TJ we do observe 
the creation of a boundary layer. The evolution can be seen in Figure [2j We remark 
that the colour range in Figure [2] is from red for 1 to blue for —1, as in Figure [TJ even 
though the extremal values for $ n during the evolution are approximately 1.03 and —1.16, 
respectively. Finally, we recall that if we choose the shape functions ( 2.11[ )(ii) or (2.11 )(iii) 
instead, then the conditions (3.11) and (3.28) yield that $ n = 1 and W n = ud for all 



n > 1 for the two schemes (3.6a b) and (3.21a b), respectively. 
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Figure 2: (e' 1 = 16 ir, ANlS°- 01) , (|2.1lfc(i 



a 



1, P = io- 



— 42, 



Creation of a boundary layer for the scheme (3.21a b). Snapshots of the solution at times 
t = 0, 8 x 10" 4 , 10" 3 , 1.5 x ID" 3 , 2.5 x 10~ 3 . 



l i- 

2' 2' 




Snapshots of the solution at times t — 0, 1, 3, 5, 7. [This computation took 6 days.] 



In the following experiments, we return to the initial data described previously, so 
that (po models a circular interface of radius Rq = 0.1. We also set p = 0. For convenience 



we recall the simulation from Barrett et al. (2012c, Fig. 5), so that (2.11) (i) applies. Here 
H = 8 and = p = 0. Moreover, = —2 and a = 0.03; and we observe that for this 
choice of parameters the condition (3.11) is satisfied if we choose e~ x = 32 7r > ^-n. A 
run for (3.6a b), with the discretization parameters Nf 
T 



4096, N c = 128, t 



32 7T > y n. 

10" 4 and 



7 is shown in Figure pj Now the advantage of the shape function choices (2.11 )(ii 



and (2.11)(iii) over the simple choice (2.11 ) (i) , as highlighted in Remark 3.3, is that 



for the same physical parameters a larger value of e can be chosen. To illustrate this, we 
repeat the simulation from Figure [3] but now for the choices (2.11) (ii) and ( ]2.11[ )(iii). This 
means that we can use e.g. e^ 1 =8ir together with the coarser discretization parameters 
Nf = 1024, N c = 128 and r = 10 -3 . The new results are shown in Figures [4] and [HJ where 
we observe the good qualitative agreement with Figure [3] We draw particular attention 
to the dramatic reduction in CPU time necessary to compute the respective simulations. 
While a further reduction in e~ l and in the discretization parameters Nf, N c and r _1 
leads to even bigger gains in computation times, the larger values of e soon lead to a loss 
of accuracy with respect to the approximation of the underlying sharp interface problem 
(1.4a-e). We illustrate this with an example for e' 1 = 2% for the choice (2.11 ) (ii) together 
with Nf = 256, N c = 64, r = 10~ 2 and T = 6. The results are shown in Figure [6j where 
we observe a qualitative difference to the three previous simulations. 



For the remainder of the simulations in this subsection we continue to employ (2.11 ) (ii) , 
but we now choose p = 0.01 A simulation corresponding to Figure [6] can be seen in 
Figure [7J We observe that in this example, the presence of kinetic undercooling (p > 0) 
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Figure 4: (e' 1 = 8vr, ANlf' 3) , (]21l|)(ii), a = 0.03, p = 0, u D = -2, ft = (-8,8 
Snapshots of the solution at times t — 0, 1, 3, 5, 7. [This computation took 4.5 hours. 




Figure 5: (e^ 1 = 8vr, ANlf' 3) , |2.1l|(iii 



Snapshots of the solution at times t = 0, 1 



a = 0.03, p = 0, m d = -2, ft = (-8,8) 2 
3, 5, 7. [This computation took 10.5 hours.] 



only has a small influence on the overall evolution. 



The remaining computations in this subsection are for the rotated hexagonal anisotropy 
ANlg. The first simulation is analogous to Figure [7j but now on the larger domain 
ft = (—16, 16) 2 . In particular, we keep all the parameters as before, apart from 7 and 
apart from Nf = 512, N c = 128 due to the increased value of H. The results are shown 
in Figure [8j We have seen in previous simulations that the value of e can have a large 
influence on the evolution of the phase field approximation. Reassuringly, in this example 
the evolution remains qualitatively unchanged if we repeat the simulation for e~ l = Air. 
A run with Nf = 1024, N c = 256, r = 10~ 3 and T = 8 is shown in Figure [9j We end this 
subsection with a repeat of the last computation, but now for the stronger supercooling 
ud = —4. The evolution now exhibits six distinct side arms, as can be seen in Figure 10 




Figure 6: {e~ l = 2vr, ANlf 3) , (|2.1l[)(ii 



a 



0.03, 



P 



0, u D 



-2, n 



Snapshots of the solution at times t = 0, 1, 3, 5, 6. [This computation took 3 minutes.] 
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Figure 7: {e~ l = 2vr, ani!°- 3) , (|2.1l[)(ii 



Snapshots of the solution at times t = 0. 



a = 0.03, p = 0.01, u D = -2, tt = (-8,8)' 
1, 3, 5, 6. [This computation took 4 minutes.] 




Figure 8: (e' 1 = 2vr, ami*, (gTTTj) (ii) , a = 0.03, p = 0.01, u D = -2, Q = (-16, 16) 2 ) 
Snapshots of the solution at times t = 0, 1, 5, 6, 8. [This computation took 25 minutes.] 




Figure 9: (e' 1 = 4vr, ami*, (2.11 ) (ii 



Snapshots of the solution at times t 
minutes.l 



a = 0.03, p = 0.01, it d = -2, tt = (-16, 16) 2 ) 
0, 1, 5, 6, 8. [This computation took 5 hours, 17 
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Figure 10: (e' 1 = 4vr, AMI*, (2.11 )(ii 



Snapshots of the solution at times t 
minutes.l 



a = 0.03, p = 0.01, u D = -4, O = (—16, 16) 2 ) 
0, 1, 2, 3, 4. [This computation took 3 hours, 5 
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Figure 11: (e' 1 = 4tt, ANlf- 3) , (|2.1l|) (ii 



Snapshots of the solution at times t 
minutes. 1 



a = 5x HT 4 , p = 0.01, u D = -0.5, Q = (-8, 8) 2 ) 
= 0, 0.1, 0.2, 0.5, 1. [This computation took 76 




Figure 12: (e 



-i 



47T, ANI 



(0.3) 



(EmRiii) 



cv 



5 x 10" 4 , p = 0.01, m d 



-0.5, ft 



(— 8,8) ) Snapshots of the solution at times t 
took 165 minutes.l 



0, 0.1, 0.2, 



0.5, 1. 



[This computation 



5.2 Stefan problem in two space dimensions 



In a first simulation for the full Stefan problem, i.e. with $ > 0, we take parameters that 
are close to the ones used in Barrett et al. (2010b, Fig. 10). In particular, we have t? = 1, 
a = 5 x 10~ 4 , p = 0.01, ud = —0.5 and R = 0.2, H = 8. An experiment with e -1 = Air 
together with Nf = 512, N c = 64, r = 10~ 3 and T = 1 is shown in Figure 11, where we 



employ (2.11 ) (ii) . We observe a very large interfacial region, which indicates that e was 



not chosen small enough. A similar behaviour can be observed for the choice (2.11)(iii), 
see Figure 12 Here we note that in this example, in line with the analysis in (2.16), 



there appears to be no benefit in using ( 2.11[ )(iii) over (2.11 ) (ii) . On reducing the size 
of the interfacial parameter e, the phase field again assumes its expected profile across 
the interface, and we obtain the following numerical results. If e~ x = 16 7r together with 

4 and T = 2 we obtain the results shown in Figure 13 
4096, N r = 128, r = 10" 4 and T 



N 



if = 2048, N c = 128, r = 10 
e^ 1 = 32 ii together with Nf = 



If 



2 we obtain the results 



shown in Figure 14 We can see that the small oscillations present in the final snapshot 



in Figure 13 have vanished in the corresponding plot in Figure 14 A closer comparison 
of the two solutions at time t = 2 can be seen in Figure 15 



The next simulation is for the rotated hexagonal anisotropy ANI3. 

1 = 16 7r together with Nf = 



parameters are as in Figure 13 
~ 4 and T 



i.e. e 



All the remaining 
2048, N r = 128, 



r = 10~ 4 and T = 2. See Figure 16 for the numerical results. The large mushy regions 
in the final plot in Figure 16 indicate once again that e needs to be chosen smaller. 
Hence we repeat this experiment and now choose e~ x = 32 n together with Nf = 4096, 
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16 7T, ANiS a3) , (|2.1l[)(ii) 



Figure 13: (e 

(— 8,8) 2 ) Snapshots of the solution at times t 
15.5 hours. 1 



a = 5 x 1(T 4 , p = 0.01, u D = -0.5, = 
0, 0.5, 1, 1.5, 2. [This computation took 




32 7T, ANlS a3) , (|2.1l[)(ii 



Figure 14: (e 

(— 8,8) 2 ) Snapshots of the solution at times t 
19.5 hours. 1 



a = 5 x 10~ 4 , p = 0.01, m d = -0.5, Q = 
0, 0.5, 1, 1.5, 2. [This computation took 





Figure 15: A close comparison of the final solutions from Figures 13 and 14 On the left 
e^ 1 = 16 Ti, and on the right e~ l = 32 tx. 
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Figure 16: (a' 1 = 16 vr, ani*, (|2.11[)(ii), a = 5 x 10~ 4 , p = O.Ol, u D = -0.5, Q = (-8, 8) 2 ) 
Snapshots of the solution at times t = 0, 0.5, 1, 1.5, 2. [This computation took 51 hours.] 










MM 



Figure 17: (e' 1 = 32 vr, ami*, (2.11 )(ii 



a 



5x 10- 4 , p = 0.01, m d = -0.5, ft 



Snapshots of the solution at times t = 0, 0.5, 1, 1.5, 2. [This computation took 19 hours.] 



N c = 128, r = 10~ 4 and T = 2. See Figure 17 for the numerical results. We observe that 
the interfacial region is now well defined and that the evolution exhibits six distinct side 
arms. 



5.3 Mullins— Sekerka in three space dimensions 



In this subsection we always employ (2.11 )(ii), and we always let $ = 0. At first we 
also choose p = 0, so that we approximate a Mullins-Sekerka problem without kinetic 
undercooling. A simulation for the cubic anisotropy ANlg with e~ l = 2tt, and for the 
physical parameters a = 0.03 and «o = —2, can be seen in Figure 18 The discretization 
parameters are Nf = 256, N c = 64, r = 10~ 3 and T = 1. We observe that the cubic 
anisotropy induces the growth of the typical six symmetric side arms. If we repeat the 
simulation with p = 0.01, which models the presence of kinetic undercooling, the shape 
of the phase field approximation of the growing crystal changes significantly. We present 
a run for the discretization parameters Nf = 256, N c = 64, r = 10 -2 and T = 2 in 
Figure [19} Note that the larger time step size used here yields a large reduction in the 



overall CPU time. 



A repeat of the simulation in Figure [19j but now for the rotated hexagonal anisotropy 

In this simulation we can observe facet breaking, both 



ani^ can be seen in Figure 20 



in the basal and in the prismal directions, similarly to the sharp interface computation 
shown in Barrett et al. (2012b, Fig. 18). With the next simulation we wish to highlight 



the effect that the choice of the mobility coefficient (3 can have on the evolution. If we 
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Figure 18: (e' 1 = 2tt, ani 9 , gTTJ(ii), a = 0.03, p = 0, u D = -2, = (-8,8) 3 ) Snapshots 
of the solution at times t = 0, 0.1, 0.2, 0.5, 1. [This computation took 13 days.] 




Figure 19: (e' 1 = 2vr, ani 9 , fl2.1l| )(ii), a = 0.03, p = 0.01, u D = -2, Q = (-8,8) 3 
Snapshots of the solution at times t = 0, 0.1, 0.5, 1, 2. [This computation took 4 days.] 



replace = 7 with /3 = /3fl a t,3, where 



is defined as in |Barrett et aT\ ( |2012b| Eq. (16)), and if we keep all of the remaining param- 
eters as before, then we obtain the results shown in Figure 21 Clearly, the growing crystal 
now assumes the shape of a flat prism. Similarly, if we choose the mobility coefficient 
= Aaii,2, where 



is defined as in Barrett et al. (2012b, Eq. (17)), then we obtain the simulation presented 



in Figure 22 This time the initially spherical crystal grows into a tall hexagonal prism. 



It is discussed in Libbrecht (2005) that different mobility coefficients are responsible for 
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Figure 20: (e' 1 = 2vr, AMI*, fl2.1l| )(ii), a = 0.03, p = 0.01, w D = -2, J] = (-8,8) 3 ; 
Snapshots of the solution at times t = 0, 0.1, 0.5, 1, 2. [This computation took 22 hours. 
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Figure 21: (e' 1 = 2tt, ami*, (2.11 )(ii), a = 0.03, p = 0.01, u D = -2, 



Snapshots of the solution at times t = 0, 0.1, 0.5, 1, 2. [This computation took 2 days.] 



the various snow crystal shapes seen in nature. In this context we remark that (1.4a 



also appears in solidification from a supersaturated solution. In this case — u is a suitably 



scaled concentration with —ud being the scaled supersaturation, see e.g. Barrett et al 



(2012b) for more details. 



5.4 Stefan problem in three space dimensions 

In this subsection we present a simulation for the full Stefan problem in three space 
dimensions for the anisotropy ANlg. To this end, we consider the physical parameters 
i? = 1, at = 10 -3 , p = 0.01, ud = —0.5 and let Q = (— 4, 4) 3 . A numerical computation 
for e^ 1 = 16 7r, together with Nf = 1024, N c = 64, r = 10 -4 and T = 0.4 can be seen 



in Figure 23 Similarly to the results in Figure 19 we observe that the growing crystal 
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Figure 22: (e" 1 = 2tt, ami*, fl2~TTj)(ii), a = 0.03, p = 0.01, u D = -2, Q = (-8,8) 3 ) 
Snapshots of the solution at times t — 0, 0.1, 0.5, 1, 1.8. [This computation took 26 
hours.] 



exhibits the typical six symmetric side arms that are common in simulations of dendritic 
growth. 
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